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ABSTRACT 

An analytical model is presented for calculating the surface density as a function of radius E(r) in 
protoplanetary disks in which a planet has opened a gap. This model is also applicable to circumbinary 
disks with extreme binary mass ratios. The gap profile can be solved for algebraically, without 
performing any numerical integrals. In contrast with previous one-dimensional gap models, this model 
correctly predicts that low-mass (sub-Jupiter) planets can open gaps in sufficiently low-viscosity disks, 
and it correctly recovers the power-law dependence of gap depth on planet-to-star mass ratio q, disk 
aspect ratio h/r, and dimensionless viscosity a found in previous numerical studies. Analytical gap 
profiles are compared with numerical calculations over a range of parameter space in q, h/r, and a, 
demonstrating accurate reproduction of the “partial gap” regime, and general agreement over a wide 
range of parameter space. 

Subject headings: hydrodynamics — planet-disk interactions — planets and satellites: formation — 
protoplanetary disks 


1. INTRODUCTION 

Disk-satellite interactions have been studied for 
decades. The general problem of an orbiting point mass 
interacting gravitationally with a surrounding disk of ei¬ 
ther gas or solids is a fundamental physics problem with a 
multitude of applications in astrophysics. Saturn’s rings, 
for example, are an exquisite testbed for disk-satellite in- 
teraction s in the case where the disk is composed entirely 
of solids (jColdreich fc TremaindllQT^ . For the hydrody¬ 
namic case, how a planet interacts with the gaseous disk 
which spawned it is of great importance to the question 
of how the planets formed. It also provides a vital key¬ 
stone in our interpretation of observational data, both of 
extrasolar planetary systems and of debris disks around 
newly formed stars. Additionally, disk-satellite interac¬ 
tions are of importance to the orbital evolution of a black 
hole binary embedded in a circumbinary disk. 

The linear theory of hydrodynamical disk- 
planet interactions has been ver y successful so far 
((Goldreich fc Tremainel [19781 II980l: IWardI 119861) . Pre¬ 
dictions of the disk’s linear response to the planet 
and subsequent predictions of the disk perturbation’s 
gravitational influence on the planet have been ac¬ 
curately reproduced in numerical studies. This is 
true for que stions of pla n etary ( Type I) migration 
rates (e.g. iTanaka et al.l I2002t iD’Angelo fc Lubo^ 
[Mnli ■ and of the eccent r icity e volution of the planet 
le.g. ICresswell & NelsonI 120061 I Cress well et al.l 120071 
i Bitsch fc S^nio lL - - 

The details of nonlinear theory, however, have not all 
been reproduced in numerical studies. In linear theory, 
a spiral wave is excited by the planet, and passes po¬ 
litely through the disk. After each wave front passes, it 
leaves no trace of its presence. A nonlinear wave, how¬ 
ever, can shock. When it does so, it imparts its angular 
momentum to the disk, causing fluid elements to move 
away from the planet’s position. This process carves out 
a low-density annulus in the vicinity of the planet’s or- 
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bit. This low-density region is known as a “gap”, and its 
detailed properties have eluded a robust and numerically 
reproducible semi-analytical description for some time. 

Many attempts have been made to describe this three- 
dimensional system with a one-dimensional model, to 
predict the gap structure, i.e. the surface dens ity E 
as a function of radius r. iVarniere et al.l (I 2004D took 
the d isk evolution equations of iLvnden-Bell fc Pringi^ 
(flTTl and added a planetary torque term, to predict 
the width and shape of gaps. The calculated gaps 
turned out to be much wide r than those seen in nu¬ 
merical studies. iCrida et al.l (j2006[l attempted to fix 
this by adding a component due to “pressure torque” in 
the disk. However, both of these studies predicted gaps 
which were much too deep. In part, this is because the 
model for torque deposition was too simplistic; for ex¬ 
ample, the torque from the planet was calculated locally, 
even though the transfer of angular momentum from the 
planet to the disk by shocks is an inherently nonlocal 
process, as the wave must propagate some distance be¬ 
fore shocking. Additionally, these studies considered the 
role of planetary torque as a barrier holding back the 
viscous accretion of gas. Suc h considerations predict 
strict criteria for gap opening (iLin fc Panaloizoul 119791 
119931: iBrvden et al.l Il999t ICrida et al.l l2006l l but semi- 
analytical and numerical studies have demonstrated that 
planets can viol ate these criteria and still open gaps on 
long t im escales (iG oodman fc Rafi kovI 120011: iDong et aP 
1201160 : iDuffell fc MacFadvenll20l3b 

Similarly problematic assumptions have been made 
for ID models of cir cumbinary disks. Disk mode ls of 
iSver fc Clarkel l)1995h and later iKocsis et al.l l)2012ll im¬ 
plicitly assume that if the planet does not migrate with 
the gas, there is a “pile-up” of gas at the outer edge of the 
gap. I n contrast, the disk model of iLubow fc D’Angelol 
(I2006H considered gas flow across the gap, but even that 
study assumed that the planet’s presence affected the 
accretion rate at infinity. 

Not only are many of these a ssumptions incorrect (see 
for example IDuffell et al.l I2014H , but they also result in 
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gap profiles which are inconsistent with numerical exper¬ 
iments. The gap profiles in these studies generally pre¬ 
dict exponentially deep gaps; that is, the surface density 
near the planet is depleted by a factor ^ exp(—(g/ggap)^) 
relative to the unperturbed state, where q is the planet- 
to-star mass ratio, and ggap is some critical threshold for 
gap opening. On the other hand, several recent numeri¬ 
cal parameter surveys have shown that gap depth scales 
as a power-law whh q, h/r and a (|Duffell fc MacFadvenI 
120131 : iFiing et al.ll201^ . so that gaps are not as depleted 
a s in the ID model s. 

iFung et m (j2014f ) provided an argument for the power- 
law scalings found in lDuffcll fc MacFadvenI l)2013f) . point¬ 
ing out that the torque produced by the planet would 
be proportional to the surface density in the gap near 
the planet’s position Ep, assuming that the torque was 
domin ated by resonances excit ed inside the gap. Very re¬ 
cently, iKanagawa et al.l (l2015| l produ ced a ID model fo r 
gaps which repeated the argument s ofiFung et al.l (120141). 
recove ring the power-law scaling of iDuffell fc MacFadvenI 

However, IKanagawa et al.l (|2015[ 1 used an incomplete 
description of the planetary torque. For most of their 
study, the torque was specified as the excitation torque 
in the disk, as calculated by linear theory, whereas 
the torque which controls gap opening is the deposition 
torque, i.e. the rate of transfer of angular momentum 
from the planetary wake to the disk by shocks. They 
attempted to fold this into their study by including a pa¬ 
rameter Xd which measured the distance the wave travels 
from the planet before shocking. The resultant gap pro¬ 
file was then some function of the quantity Xd, which was 
treated as a free parameter. 

Fortunately, the process of excitation and depo- 
sition of torque in the disk has been studied by 
iGoodman fc RafikovI (j2001[) (hereafter GROl). There¬ 
fore, there is no need to speculate about this pro¬ 
cess, at least in the weakly nonlinear regime. In the 
present study, the planetary angular momentum flux de- 
rived by GROl i s app lied to a disk model similar to 
IKanagawa et al.l (120151 1. These disk equations can then 
be solved algebraically for the surface density, so that it 
is not necessary to perf orm any numerical integ rals as in 
iCrida et"^ (1200611 and IKanagawa et al.l (l2015ll . 

The model is described in section [2j In section [3l the 
resultant gap profile E(r) is co mpared with gap profiles 
from a recent parameter survey (jDuffellll2014[l and agree¬ 
ment is found, particularly in the “partial gap” regime. 
This provides confidence that a complete ID disk model 
may be possible after taking into account higher-order 
effects. A summary and discussion is given in section |4l 

These results may be testable observationally in the 
not-too-distant future, as the ALMA project may be ca¬ 
pable of measuring surface density profiles in transition 
disks. 

2. DESCRIPTION OF THE ANALYTICAL MODEL 

The i nitial approach here i s essentially equivalent to 
that of iVarniere et al.l (|2004f l. though the underlying 
equations will be derived in a slightly different manner. 
Rather than using the equation of mass flux to derive a 
differential equation for S(r), a planetary torque term 
is inserted into the equation of angular momentum flux, 
and the fluxes at infinity are assumed to be unaffected 


by the presence of the planet. The system is then solved 
algebraically. 

The equations of mass and angular momentum con¬ 
servation in a steady-state viscous disk are as follows 
(ignoring the planet, at first): 

M = —2'KrTjVr = const (1) 

j = = const (2) 

dr 

where M is the mass accretion rate, J is the angular 
momentum flux, Vr is the radial drift velocity, is the 
orbital angular velocity of the gas, and v is the kinematic 
viscosity. All quantities are vertically and azimuthally 
averaged. The accretion rate is defined so that positive 
M implies inward accretion. The first term in ([2]) is the 
flux due to advection, and the second term is due to 
viscous torques. 

Assuming Keplerian D (as will generally be done for 
simplicity in this study), equation ([2|) becomes 

j = —Mr^VL -I- 37rE^r^O. (3) 

There are two possible solutions such that M and J 
are both spatially uniform. The first is M = 0, J = 
iTrYjvr^Q. = constant, and the second is M = SttS:/ = 
constant, J = 0, so that viscous angular momentum 
transport outwards exactly cancels advective transport 
inwards. The latter solution will be concentrated on in 
this study, as there is zero accretion in the former solu¬ 
tion. However, the M = 0 solution does not pose any 
additional challenges, and what follows would be nearly 
identical (a linear combination of the two is also possi¬ 
ble, but again this is ignored in the current study for 
simplicity). 

Thus, for this (unperturbed) solution, the surface den¬ 
sity E is inversely proportional to the viscosity, v. If ^ is 
spatially uniform, then E will also be, in steady state. 

Once the constants M and J are known, one can add 
a planetary torque as an additional source of angular 
momentum flux: 

j = —Mr^{l -I- 37rE^r^O -|- ‘&p(r) = 0 (4) 

<l>p(r) is the angular momentum flux generated by the 
planet; the torque density deposited by shocks is given 
by the negative derivative of this function. 

Note that it is assumed that the planet does not affect 
the values of M and J at infinity, so that in particular 
the gap does not act as a “barrier” preventing mass flux 
through the system. The surface density E(r) adjusts 
itself to maintain these values of M and J in the presence 
of the planet. The value of M is the unperturbed value, 
the same as it would be in absence of the planet: 

M = ZirvYi^ir) (5) 

(the viscosity v can also depend on radius, but this no¬ 
tation was omitted for simplicity). The planet’s influence 
is entirely housed in the planetary angular momentum 
flux, 4>p(r). 4>p(r) is determined by how far the wave 
propagates before shocking, and how quickly its angu¬ 
lar momentum is transferred to the disk after shocking. 
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Fig. 1.— Gap depth Sp/Sp is plott ed for 84 differ ent planet-disk 
models, from the parameter survey of lDuffelll 1120141 1. These depths 
are compared with equation (9). Planet-to-star mass ratio ranged 
from q = 10“® to q = 2x 10“®, Mach number ranged from At = 14 
to At = 40, and viscosity ranged from a = 0.005 to a = 0.04. The 
horizontal axis plots the quantity if ( 9 ) = ja. 


This p rocess has been studied by iGoodman fc RafikovI 
(|2001ll , and their calculations provide a clear prescription 
for ‘hp('r) which applies in the weakly nonlinear regime. 
The specific details of this function will be discussed fur¬ 
ther in section [2.11 but the overall scaling is given by the 
one-sided torque emanating from the planet: 

<l>p(r) = (6) 

where a is the radial position of the planet, Ep is the 
surface density at the planet’s position E(a), q is the 
planet-to-star mass ratio, and Ai is the Mach number 
(the inverse of the disk aspect ratio h/r). f{r) is a dimen¬ 
sionless function of r, discussed later. Near the planet, 
/ evaluates to /(a) = /o, where /o is some dimensionless 
order-unity constant. Note that there is an implicit as¬ 
sumption that the waves are excited within the gap, so 
that the angular momentum flux is proportional to the 
density in the gap, Ep. 

Evaluating dl} at the planetary position (also incorpo¬ 
rating m for the mass flux) gives: 


— STTz/Ega^flp -I- STrEpZ/a^rip -|- fo = 0. (7) 

This can be inverted to find a formula for gap depth: 


Ep fT , j. <l^A4^a^flp I 
^ = + Sttzz ) 


( 8 ) 


This recovers the gap d epth scaling reported in 
iDuffell fc MacFadvenI (12011^. This exercise was also 
effectively carried out by iFung et al.l l)2014ll . and by 
iKanagawa et al.l l)2015f) . More specifically. 


i_ (9) 

Eo 1 + foK{q)/{3TT) ^ ’ 

where 

K{q) = q'^M^/a (10) 

and a is the dimensionless Shakura-Sunyaev viscosity 
parameter, a = j Figure [T] compares gap 
depths from equation (jl|) with steady-state 2D solutions 
from a nume rical parameter survey over q, a and M 
(|Duffellir2014ll . demonstrating the accuracy of this for¬ 
mula over a wide range of parameter space. For these 
84 disk models, the only significant deviation from equa¬ 
tion ([5]) is in the K ^ 10 — 1000 range, where higher- 
or der effects become important. For larger values of 
A", iFung et'all (|2014f) found that this scaling eventually 
changes, suggesting that higher-order effects become im¬ 
portant for very large K > 10^. This figure was used to 
constrain the one free parameter /o = 0.45. This value of 
/o is also consistent with the one-sided torque calculated 
numerically in these disk models (/o « 0.451 ±0.534/Al, 
where the ± sign indicates the torque excited on either 
side of the planet, with the greater torque excited in the 
outer disk). After choosing /o = 0.45, the model is com¬ 
pletely specified. 

2.1. Density as a Function of Radius 

Now with this value of Ep, one can return to the orig¬ 
inal J equation: 


— 37ri/Eo(r)r^O -|- 37rE(r)izr^D -|- Jlpa'^D.y f (r) = 0. 

( 11 ) 

Given T,p/llo, which was calculated above (O, one can 
solve (fTTll for E(r): 


E(r) 


Eo(r) 


f{r)K{q)/{3TT) 

1 + /oA:(g)/(37r) 



( 12 ) 


The function /(r) is a scaled-out version of the an¬ 
gular momentum flux due to the shocking of the plan- 
eta ry wake, which has b e en cal culated semi-analytically 
by IGoodman fc RafikovI (j2001l l in shearing-box coordi¬ 
nates, and thi s result has be en extended to global disk 
coordinates bv IRafiko^ (|2002ll . Therefore, there is a well- 
specified formula for 4>p(r) throughout the entire disk, 
specifically applicable for low-mass planets whose plan¬ 
etary influence is weakly nonlinear. By design, this is 
only meant to apply if the quantity qA4^ <C 1, meaning 
for sub-Jupiter planets. However, it may be possible to 
improve the semi-analytic calculation of ‘f’p(r), leading 
to an accurate model for all planet masses. This will be 
attempted in a future study. 

For the present stud y, the formula for <h „ (r) c omes 
from the calcu lations of lGoodman fc RafikovI (1200111 and 
I RafikovI (|2002ll . The wave conserves its angular momen¬ 
tum until it shocks and deposits its angular momentum 
into the disk. Therefore, the function f(r) is simply a 
constant, until the distance from the planet that the wave 
shocks. The explicit shape of this function is shown in 
Figure 3 of GROl, and its extension to a global disk is 
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q = 3.2x10'^ q = 6.4x10''’ 








Fig. 2.— Several planet masses test the ID model, for a disk 
with A4 = 20, a = 0.01. Solid curves are numerical, while dashed 
curves are analytical. Agreement is found for lower-mass planets, 
in the “partial gap” regime. For more massive planets {q ~ 10~^, 
or Jupiter-mass), there is less agreement. However, improvements 
may be made by taking into account the change in the torque due 
to the detailed gap structure. 

shown in Figure 3 of iRafiko^ (|2002f) . The function /(r) 
can be described in terms of the parameter T(r): 


fir) 


fo _ T{r) < Tsh 

fo\/Tsh/T{r) r(r) > r^h 


(13) 


where the shock positi on Tsh was calculated by 
iGoodman fc R.afikovI (j2001h . and is given by 


= 1.89 + 0.53/(gA43) (14) 

The parameter r(r) represents an appropriately-scaled 
distance from the planet. In shearing-box coordinates, 
this is given by 


p3/4 q 

r{r) = —\-M{r/a-l)f/^ (15) 

To apply this to a global disk, in which surface density 
and sound speed vary as S oc r~P^ and c oc one 

must perform the following integral 



Fig. 3.— Variation with viscosity is tested for the g = 6.4 X 10 ® 
(~ 2 OM 0 ) planet, in a disk with A4 = 20. 


T(r) 




1^3/2 _ i|3/2^(5p.+pr)/2-ll/4^^ 


which can be integrated analytically and expressed in 
terms of hypergeometric functions. 

Equations (USD, (HD), and (HSl) completely specify 
the value of the surface density at every point in the disk, 
given the parameters g, a, and h/r, and the background 
surface density profile Eo(r). 


3. COMPARISON WITH NUMERICAL RESULTS 

The accuracy and range of applicability of this simple 
analytical model are tested by comparing directly with 
numer ical calculation s from the extensive parameter sur¬ 
vey of iDuffell (j2014f ). In that survey, three parameters 
were varied: mass ratio g, Mach number M (or equiv¬ 
alently h/r = 1/A4), and viscosity v (or equivalently, 
a = / {afVLp). 

Figure [2] shows variation with mass ratio for a disk 
with A4 = 20 and a = 0.01. Mass ratios from q = 
3.2 X 10“^ to q = 10“^ are studied showing a wide range 
of applicability of the ID model. Near q = 10“^ (Jupiter- 
mass), the model begins to break down as the gap is 
significantly wider than predicted, but the depth of the 
gap is still reasonably accurate. 

The reason for the deviation around Jupiter-mass is 
likely that all of the excitation is assumed to come from 
within the gap, whereas once the gap becomes deep 
enough, the dominant resonances live just outside the 
gap edges. It is likely that a great deal of improvement 
would come from modifying the function 4>p(r) to ac¬ 
count for excitation of waves in the presence of the de¬ 
tailed shape of the gap. 

Figure [3] shows variation of the gap profile with vis¬ 
cosity, for a g = 6.4 x 10“® ('^ 2 OM 0 ) planet in a disk 
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h/r = 0.071 h/r = 0.05 






Fig. 4.— Dependence on disk temperature is explored, by vary¬ 
ing the Mach number J\4 (or equivalently hjr = l/Ai). a. is kept 
fixed, whereas v = is varied. 

with h/r = 0.05 {M = 20). The gap depth and shape is 
accurately captured for this range of viscosities. 

Finally, Figure |4] shows dependence on disk tempera¬ 
ture (equivalently M or h/r) for g = 6.4x 10“®, a = 0.01 
{a is fixed, but v = aa?^p/M^ varies with M). The 
shape of the gap is reasonably well-captured, but again 
once the gap becomes deep enough, there is a similar 
deviation from the model (this same deviation occurs at 
around AT ^ 100 in Figure [1] note that the last panel in 
this figure has K = 42). 

This analytic model should be improved for the case 
of deep gaps (this is planned for a future study), but 
it accurately captures the “partial gap” regime, and the 


basic fundamentals a r e corr ect in that the scalings of 
iDuffell fc MacFadvenI (|2013f) are reproduced, and low- 
mass planets can also open gaps in sufficiently low- 
viscosity disks, wh ic h has also been seen num erically 
(jPong et al.ll2011bl l^: IDuffell fc MacFadvenll201^ . 

4. DISCUSSION 

A simple analytical model has been presented for gap 
profiles in protoplanetary disks. In this model, effects 
due to planetary torque come directly from the detailed 
calcul ations of shock propagation by GROl, and iR.afikovI 

IMl). 

The model is solved algebraically, and compared with 
long-duration numerical calculations. The model’s per¬ 
formance in the “partial gap” regime is unprecedented in 
its agreement with numerical calculations and its repro¬ 
duction of empirically found scalings. For deep gaps, the 
model does not correctly predict the gap shape. This is 
most likely due to the fact that resonances located several 
scale heights from the planet become dominant once the 
gap suppresses the resonances nearest the planet. Ad¬ 
ditionally, effects due to non-Keplerian orbital motion 
should become important for very deep gaps. 

One could envision an iterative procedure for folding 
these effects into the analytic formula. For example, 
given the profile S(r), one could calculate the contribu¬ 
tion of the torque from all resonances in the disk. This 
could give a new form for $p(r’), which would provide a 
corrected surface density via equation (ED- This process 
could be repeated until a converged gap profile is found. 
Such a procedure will be attempted in a future study. 

Resources supporting this work were provided by the 
NASA High-End Computing (HEC) Program through 
the NASA Advanced Supercomputing (NAS) Division at 
Ames Research Center. 

I am grateful to Eugene Chiang, Roman Rafikov, 
Andrew MacEadyen, Steve Stabler, Eliot Quataert, 
Kazuhiro Kanagawa, Alessandro Morbidelli, and Aure- 
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